%Script to produce visualizations of DOW radar reflectivity data during the
%Marshall Fire

%Produced by: Neil P. Lareau summer 2022; for questions contact
%nlareau@unr.edu

close all;
clear all;

%% Load the data (this comes from the associated script: "marshall_radar_volumes.m")
%you can load either "marshall_plume_d1.mat" for the first deployment site
%or "marshall_plume_d2.mat" for the second deployment site

load('marhsall_plume_d1.mat','dbzmaster','timemaster','xxx','yyy','zzz','stationlat','stationlon','stationalt');
%these data are the 3D cartesian gridded post-processed version of the
%"raw" DOW observations. The 3D grids allow the production of volume
%renderings.

% set up grid to house the data
xv=(-2*10^4):100:(4*10^4);
yv=(-1.6*10^4):100:(.5*10^4);
zv=[0:100:5000];

%% Load the Marshall Fire Perimeter JSON file 'marhsall.json'
fileName='marshall.json';
str = fileread(fileName); % dedicated for reading files as text
test=jsondecode(str);

%produce a visualization to make sure it gives the expected result.
figure(1);clf;
for aa=1:size(test.features.geometry.coordinates)
    figure(1);hold on;
    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot(pnow(:,1),pnow(:,2),'r')
end

%% Load Terrain Data using tiles of SRTM 30 meter data
lons=stationlon;
lats=stationlat;
re=6.37*10^6; %earth radius
scalefactor=(re.*cosd(nanmean(lats))).*(pi./180); %scale factor to map dlon to dx

addpath('~/Dropbox/MATLAB/')
X1=readhgt('~/Dropbox/MARSHALL_FIRE/N40W106.hgt');
[ltmt1,lnmt1]=meshgrid(X1.lat,X1.lon);
ltall=ltmt1(1:5:end,1:5:end);%
lnall=lnmt1(1:5:end,1:5:end);%
Zall=X1.z(1:5:end,1:5:end)';%;
terrain_lon=lnall;%mt1;
terrain_lat=ltall;%mt1;
xterrain=double((terrain_lon-mean(lons)).*scalefactor);
yterrain=double((terrain_lat-mean(lats)).*111180);
Zterrain=double(Zall);

%
X2=readhgt('~/Dropbox/MARSHALL_FIRE/N39W106.hgt');
[ltmt2,lnmt2]=meshgrid(X2.lat,X2.lon);
ltall2=ltmt2(1:5:end,1:5:end);%
lnall2=lnmt2(1:5:end,1:5:end);%
Zall2=X2.z(1:5:end,1:5:end)';%
terrain_lon2=lnall2;
terrain_lat2=ltall2;%
xterrain2=double((terrain_lon2-mean(lons)).*scalefactor);
yterrain2=double((terrain_lat2-mean(lats)).*111180);
Zterrain2=double(Zall2);

X3=readhgt('~/Dropbox/MARSHALL_FIRE/N39W105.hgt');
[ltmt3,lnmt3]=meshgrid(X3.lat,X3.lon);
ltall3=ltmt3(1:5:end,1:5:end);%cat(1,ltmt1,ltmt2);
lnall3=lnmt3(1:5:end,1:5:end);%cat(1,lnmt1,lnmt2);
Zall3=X3.z(1:5:end,1:5:end)';%cat(1,X1.z',X2.z');
terrain_lon3=lnall3;%mt1;
terrain_lat3=ltall3;%mt1;
xterrain3=double((terrain_lon3-mean(lons)).*scalefactor);
yterrain3=double((terrain_lat3-mean(lats)).*111180);
Zterrain3=double(Zall3);

%
X4=readhgt('~/Dropbox/MARSHALL_FIRE/N40W105.hgt');
[ltmt4,lnmt4]=meshgrid(X4.lat,X4.lon);
ltall4=ltmt4(1:5:end,1:5:end);%cat(1,ltmt1,ltmt2);
lnall4=lnmt4(1:5:end,1:5:end);
Zall4=X4.z(1:5:end,1:5:end)';
terrain_lon4=lnall4;%mt1;
terrain_lat4=ltall4;%mt1;
xterrain4=double((terrain_lon4-mean(lons)).*scalefactor);
yterrain4=double((terrain_lat4-mean(lats)).*111180);
Zterrain4=double(Zall4);


%% Now produce 3D visualizations

dbzlist=[-10 10 12 17 21 26 29]; %list of dbz surface values
clrs=[.7 .7 .7; .7 .7 .7; 1 .7 .3; 1 .5 .1; 1 0 0; .9 0 0; .8 0 0]; %list of face colors for plumes
trans=[.1 .2 .3 .3 .5 .7 .7]; % surface transparency
alt=stationalt; %altitude of the DOW

%% find a target time for producing 3D visualizations of the plume
[~,tidx]=min(abs(timemaster-datenum(2021,12,30,22,22,32))); %you can adapt the time in the "datenum" call as you see fit

for vv=tidx % this was set up as a loop at first, but in its current form just uses the above defined target time
    %open figure and give it a size
    figure(200);clf; hold on;set(gcf,'color','w','pos',[100 100 1200 800])
    hold on;
    strd=1; %makes plotting faster, but can be set to 1 to get full terrain resolution if needed
    surf(xterrain(1:strd:end,1:strd:end),yterrain(1:strd:end,1:strd:end),Zterrain(1:strd:end,1:strd:end));shading flat;
    surf(xterrain2(1:strd:end,1:strd:end),yterrain2(1:strd:end,1:strd:end),Zterrain2(1:strd:end,1:strd:end));shading flat;
    surf(xterrain3(1:strd:end,1:strd:end),yterrain3(1:strd:end,1:strd:end),Zterrain3(1:strd:end,1:strd:end));shading flat;
    surf(xterrain4(1:strd:end,1:strd:end),yterrain4(1:strd:end,1:strd:end),Zterrain4(1:strd:end,1:strd:end));shading flat;

    daspect([1 1 .4]); %set the vertical stretching using 0-1 in the last place
    cmp=flipud(bone(64))*.5; %dim gray colormap for terrain
    colormap(cmp)
    caxis([800 2300]); %color limits for terrain near boulder
    hold on;
    material dull;
    grid on; box on; set(gca,'linewidth',2,'layer','top');
    camlight(40,10,'infinite')
    lighting flat;
    hold on;
    grid on; box on; set(gca,'linewidth',2,'layer','top');
    lighting gouraud;
    zlabel('Altitude [m]','fontsize',15,'fontweight','bold');
    xlabel('x-dist [m]');
    ylabel('y-dist [m]');
    set(gca,'fontsize',15,'fontweight','bold')



    %create transparent isosurfaces at the dbz list values
   
    % in this call we take the time-maximum reflectivity for every grid
    % point in a time window of +/- 5 time steps (which translates to ~10
    % minutes of elapsed time in the DOW scans
    dbznow=squeeze(nanmax(dbzmaster(vv-5:vv+4,:,:,:),[],1)); %you can change the time-maximum window as you see fit

    %% trim data to remove spurious ground returns if necessary
    dbzmap=nanmax(dbznow,[],3);%first make a map of the column maximum reflectivity
    figure(51);clf; hold on; %plot these data
    pcolor(squeeze(xxx(:,:,1)),squeeze(yyy(:,:,1)),dbzmap);shading flat;
    figure(51);
    [xe,ye]=ginput(); %draw a polygon using mouse clicks to encompass the ground returns (which are west and NW of the main plume)
    
    iin=inpolygon(squeeze(xxx(:,:,1)),squeeze(yyy(:,:,1)),xe,ye); %find grid points that are in the hand drawn polygon
    inn3=repmat(iin,1,1,size(yyy,3));%make a 3d grid of the "inpolygon" points
    dbznow(inn3==1)=-30; %set refelctivity values in this region to a clear air value of -30 dbz
     

    %% Produce the 3D isosurfaces
    figure(200)
    for dbs=1:length(dbzlist) %loop over reflectivity values
        %now use a 3D smoothing function with different spatial resolution
        %as we move toward lower reflectivity values
        if dbzlist(dbs)>=20; 
            p = patch(isosurface(xxx,yyy,zzz+nanmean(alt),smooth3(dbznow,'gaussian',[9 9 9]),dbzlist(dbs)));
        elseif dbzlist(dbs)==-10
            p = patch(isosurface(xxx,yyy,zzz+nanmean(alt),smooth3(dbznow,'gaussian',[11 11 11]),dbzlist(dbs)));
        else
            p = patch(isosurface(xxx,yyy,zzz+nanmean(alt),smooth3(dbznow,'gaussian',[15 15 15]),dbzlist(dbs)));
        end
        p.FaceColor = clrs(dbs,:); %set the face color of the surface
        p.FaceAlpha=trans(dbs); %set transparency
        p.EdgeColor = 'none'; %turn off edge colors
    end

    lighting gouraud; 
    camlight(10,-400);%define camera light position
    daspect([1 1 .4]); %vertical strecting. 
    grid on; box on;
    set(gca,'fontsize',18,'fontweight','bold','linewidth',2,'layer','top')
    xlabel('x-dist [m]');
    ylabel('y-dist [m]');
    zlabel('z-dist [m MSL]')
    xlim([-10*10^3 2.5*10^4]);%max(xv)])
    ylim([min(yv) max(yv)])
    zlim([1000 5000])
    title(strcat(datestr(timemaster(vv-5),'HH:MM'),'-',datestr(timemaster(vv+4),'HH:MM'))) %title with span of times
    view(-63,37); %set the camera view

    %% ADD the fire perimeter data along the sloping surface

    for aa=1:size(test.features.geometry.coordinates)
        %     figure(1);hold on;
        pnow=[squeeze(test.features.geometry.coordinates{aa})];
        plot3((pnow(:,1)-nanmean(lons)).*scalefactor,(pnow(:,2)-nanmean(lats)).*111180,1800*ones(size(pnow(:,2))),'r','linewidth',2)
    end

    export_fig marshall_d2a_times.tif -append -nocrop
    pause(.1)
    
end


%% Create a rotating animation of the volume rendering by moving the azimuth and elevation of the view and appending the images to a tif
for azz=214:2:500
    view(azz,1+azz/10)
    pause(.4)
    export_fig marshall_fire_rotate_AMS.tif -append -nocrop
end


%% convert to animated gif
im2gif('marshall_fire_rotate_AMS.tif','-delay',.2)

%% Generate terrain cross section
xer=[-2 4 4 -2 -2]*10^4;%max(xv)])
yer=[-8 -8 0 0 -8]*10^3;
i1=inpolygon(xterrain,yterrain,xer,yer);
i2=inpolygon(xterrain2,yterrain2,xer,yer);
i3=inpolygon(xterrain3,yterrain3,xer,yer);
i4=inpolygon(xterrain4,yterrain4,xer,yer);

xbig=[xterrain(i1); xterrain2(i2); xterrain3(i3); xterrain4(i4)];
ybig=[yterrain(i1); yterrain2(i2); yterrain3(i3); yterrain4(i4)];
zbig=[Zterrain(i1);Zterrain2(i2);Zterrain3(i3);Zterrain4(i4)];
FT=TriScatteredInterp(xbig(:),ybig(:),zbig(:));
[xx,yy]=meshgrid(xv,yv);

Tint=FT(xx,yy);
Tint=reshape(Tint,size(xx));

figure(8);clf;
subplot(3,1,1:2);cla;
pcolor(xx,yy,Tint);shading flat; caxis([0 3000])
subplot(3,1,3);cla
plot(xv,squeeze(nanmax(Tint,[],1)));

%% Create a cros section (across the wind) of the radar reflectivity
for ttt=tidx %use same time as defined above for the 3d image

    load('radar_cmap.mat') %this is a convectional radar colormap
    figure(5);clf; set(gcf,'color','w','pos',[100 100 1000 300]) %call and define figure props
    dbzbar=squeeze(nanmax(nanmax(dbzmaster((ttt-5):(ttt+4),:,:,:),[],2),[],1)); %this is the time-maximum, meriodional maximum reflectivity
    dbzbar(dbzbar<-10)=nan; %blank data less than -10 dbz
    pcolor((xv./scalefactor)+stationlon,zv+nanmean(alt),medfilt2(dbzbar,[3 3])');shading flat; %pcolor mesh with 3x3 median filter
    colormap(cmap)
    caxis([-20 30]);
    cbh=colorbar; ylabel(cbh,'dbZe','fontsize',15,'fontweight','bold')
    grid on; box on; set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
    xlabel('Longitude')
    ylabel('Alt [m MSL]')
    pause(.1)
    hold on;
    ah=area((xv./scalefactor)+stationlon,squeeze(nanmax(Tint,[],1)));%,'k','linewidth',5);
    ah.FaceColor=[.7 .7 .7];
    ah.LineWidth=3;
    ah.EdgeColor='k';
    ylim([1200 6000])
    title(strcat(datestr(timemaster(ttt-5),'HHMM'),'-',datestr(timemaster(ttt+4),'HHMM')));
%     export_fig MARSHALL_CROSS_SECTION.png -m2.5
end
